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A RAPID COMPUTATIONAL MODEL FOR ESTIMATING THE 
PERFORMANCE OF COMPLIANT AIRFOILS IN CASCADES 

F. Sisto and M. Avila 
SUMMARY 

By designing the blades in a turbomachine to have a specific 
schedule of structural stiffness (typically more compliant than 
normal) it is possible to endow the end product with an aero- 
thermodynamic parfomnance characteristic that is different from the 
rigid-bladed counterpart. Presumably the resulting characteristic 
of the compliant turbomachine, within certain operational limits, 
is preferable to the more rigid turbomachine. Since there is no 
active controller involved in changing the blade geometry in 
response to changing operating conditions, the process may be 
thought of as passive control by aeroelastic tailoring. In fact, 
since the objective function is a change in the performance 
characteristic, the overarching process may be termed '’aerothermo- 
elastic tailoring": the principles of static aeroelasticity are 
applied to a turbomachine to change its performance in a thermo- 
dynamic sense. 

The problem posed is extremely complex stemming from the need 
to couple the flow computation over a deformable surface with the 
structural dynamics computation of the body occupying the other 
side of the internal boundary. If no aeroelastic instabilities are 
encountered at a particular operating point the resulting con- 
figuration is fixed in time and the static solution (flow and 
structure) is a special case of the more general dynamic 
formulation. 

In this situation it is useful to develop an approximate model for 
use in the rapid exploration of many static aeroelastic configura- 
tions, described by a host of controlling parameters. In this 
fashion those regions of parameter space showing potential 
advantages for the compliant-bladed turbomachine may then be 
computationally refined using more accurate, but more time- 
consuming, CFD coupled to discretized structural dynamic codes 
(FEM) . It is the purpose of the present model to permit this rapid 
exploration of parameter space to delineate the more promising 
regions for the subsequent detailed studies. 

I. AERODYNAMICS 



For the rapid computation of the aerodynamic loads acting on 
airfoils in cascade, a two-dimensional incompressible potential 
flow model is chosen using a vortex sheet representation of the 
thin airfoils. Under these assumptions the upwash induced along the 
chord of the reference (zeroth) airfoil is given by 



( 1 ) 



1 1 

2 J w 2 » 

0 ^ 0 

Here the variables x(and $) are chordwise coordinates made dimen- 
sionless w.r.t. the chord, c, and g(x) is the unknown strength of 
the vortex sheet representing each airfoil. The first integral 
represents the contribution of the reference airfoil and the second 
integral represents the contribution from the remainder of the 
cascade through the "influence function" F under the assumption 
that every blade has the identical distribution g(x). 

Figure 1 shows the geometric configuration under analysis. Eq 
(1) is derived in Appendix A, as well as the inversion of the 
integral equation for the solution g(x) under appropriate non- 
penetration boundary conditions and specified vector mean velocity 
of flow (V,flt) . 

The solution is obtained under the classical thin airfoil 
coordinate transformation 



x-Vi (1-COS0) 



( 2 ) 



6-COS'^ (l-2x) (2') 

followed by a discretization of the transformed variable 

(3) 



where n is the number of segments along the chord, or n+1 is the 
number of gridpoints at which the boundary conditions are satisfied 
and also the values of g(6) are determined. 

Using results first obtained in (Sisto, 1967) the distributed 
circulation can be expressed in the following form 



g(x)-V, 






1-x 



f(x) 



(4) 



(thus rendering f dimensionless) and the solution for f in matrix 
format appears as follows: 



f (6) -u;- -Huj 



(5) 
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Here A is the (n+l)x(n+l) aerodynamic matrix and U is the upwash 
column matrix or vector. 

The expression for the upwash U depends on the mean flow angle 
a and y' the camberline slope distribution. The expression for A is 
quite complex. It reflects the geometry of the cascade: s/c, ^ and, 
again, y' as a result of the method of inverting the fundamental 
integral equation. The relevant expressions appear as Eqs. (A2) 
(A14) and (A16) determining A and Eqs. (A14) (A19) (A20) and (A21) 
for U. Both have been programmed for solution in the Matlab 
language and form part of the program testOO.m and solve. m 
appearing as Appendix C. 

II. STRUCTURE 

Each of the identical airfoil sections comprising the cascade 
are modelled as two-dimensional plates of small curvature and 
variable thickness. Each plate may be thought of as the character- 
istic or representative section of an actual blade at 75 or 80 
percent of the span of a typical cantilever blade. As such the 
elastic displacements may be expressed in terms of the bending 
moment, M according to 



d^y _ M{x) 
dx2 ' EU) 



( 6 ) 



where E is the local stiffness related to Young's modulus, cross- 
sectional moment-of-inertia and Poisson's ratio if a solid 
orthotropic material is under consideration. For an aeroelastically 
tailored structure it is presumed that E is under the control of 
the designer within limits determined by strength considerations, 
stability, etc. 

Recognizing that changes in stagger related to the untwist of 
a three-dimensional blade may be an important aeroelastic para- 
meter, the two-dimensional counterpart model will be supported at 
the leading-edge and near the trailing edge by springs whose 
stiffness may be specified. In actual construction of a compliant 
blade these springs represent the restraint to bending provided by 
the spanwise spars required to support the more compliant material 
comprising the airfoil shaped profile. 

The elastic bending moment M(x) is therefore a function of the 
aerodynamic loading whereas the elastic displacements depend as 
well upon the strength of leading edge and rearward springs. The 
structural configuration is sketched in Fig. 2 where only a single 
airfoil need be shown, the others being identical. (The case of a 
shrouded cascade of airfoils would need additional modification of 
the model) . 

The aeroelastic coupling is provided by this aerodynamic 
loading acting on the compliant structure. (There may be a small 
amount of cambor in the unstrained state upon which the elastic 
displacements are superposed) . The aerodynamic load is Jjrovided by 
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dL - CVdT - liVgix) d{cx) 



dL - CV\ f(x)d(cx) 



( 7 ) 



N ^ 



It is this loading which is applied to the plate; the details are 
developed in Appendix B. It may be noticed that dL has a square 
root singularity at the leading edge. However, under the coordinate 
transformation given by Eq (2), the singularity is integrable in 
the determination of lift and moment, viz 



is the dynamic presssure based on the vector mean velocity, V. For 
known aerodynamic loading and structural stiffness, Eq (6) may be 
solved for the elastic displacements Y(x). Subsequently the un- 
strained slopes may be augmented by the elastic slopes to obtain 
the overall slope distribution. This total slope function comprises 
the input to the aerodynamics subroutine for computation of the 
aerodynamic loading. 

III. COMPUTATIONAL LOGIC 

Conceptually the problem at hand is an eigenvalue problem 
similar in concept to the problem of wing divergence. However, 
owing to the complexity of the aerodynamic formulation, the method 
of solution is by iteration. An assumed blade shape produces a 
specific aerodynamic load. This load applied to the blade structure 
produces a specific schedule of displacements. These displacements 
define a new blade shape and completes the first loop of the 
iteration. If the iterations converge the stable static solution is 
obtained. If the iterations do not converge the condition of 
"divergence" states that the initial set of aerodynamic and 
structural parameters are in supercritical region of parameter 
space and no stable solution is obtainable. 

The looping concept is illustrated in Fig. 3, a flow chart of 
the computational logic used to iterate the solution. 

This chart, along with the actual program statements (ex- 
pressed in MatLab language) comprise Appendix C. 

IV. PARAMETRIC STUDY 

A preliminary use of the completed program demonstrated was in 
the brief parametric study described below. The program solve. m 
(with input provided from a separate file testOO.m) takes the 
converged result from the loop and computes a number of related 
parameters . 



dL - gc(l+cos0) /(6) d0 




where 



q - VzCV^ 



( 8 ) 
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In order to assess the proximity to stall a very elementary 
approach is taken. The airfoil cascade is assumed to develop the 
deflection (lift) and drag provided by Howell's crude and early 
correlation of the cascaded airfoil. Thus the camber shape, stagger 
and solidity can be used to enter Howell's charts and predict the 
profile drag and lift (deflection) coefficients. 

These data are developed for the representative section of 
single rotor compressor at the radius, or spanwise location, of the 
representative or characteristic section. For this study a design 
point wheel speed u=468 ft/sec, stagger angle 0 = 45 ° and unit 
solidity are chosen* and various graphs of loading coefficient T 
and rotor efficiency, rj, can be plotted against flow coefficient, 
0. These quantities derived from Howell's correlation are not 
presented as being quantitatively accurate. However, they do 
represent a consistent system of predicting drag-related 
performance para-meters, and the relative values are thought to 
indicate the 

proximity to stall, or extent of separation. Thus, the rotor 
efficiency at the same flow coefficient is a measure of stall 
mitigation if the efficiency of the compliant rotor is higher than 
the efficiency of the rigid rotor. 

In Fig. 4 and the charts following the results of an 
alternative study are shown for the rigid and compliant rotors 
having the fol owing characteristics 



Unstrained ''amber, deg. 
Unstrained Stagger, deg. 
Solidity 

of Elasticity 
Design wheel speed, fps 



Rigid Blades Compliant Blades 



22 

45 

1.0 

1 . SxlO^psi 
468 



15 

45 

1.0 

1 . OxlO'^psi 
468 



In this study the compliant and rigid blades have been tailored so 
as to have virtually the same performance at the wheel speed, u=4 68 
fps. There the system is assumed to have an important second 
operating condition at 10% overspeed, or U=515 fps. The performance 
at this second condition then discriminates between the two rotors 
vis a vis stall mitigation. 



Since the aerodynamics formulation is inviscid in the interests of 
rapidity of computation, an alternative assessment of the propen- 
sity to stall is made on the basis of several measures: incidence 
and a leading edge loading parameter taken in conjunction with the 
camber angle. The last is expected to change significantly with 
changing operating conditions for the compliant blade. 
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(lower) for the compliant rotor at all flow coefficients in the 
overspeed condition. It should be recalled that both rigid and 
compliant rotors have been configured to have virtually identical 
geometry and performance at design speed. 

The difference shown in Fig. 4 has not been optimized, or 
maximized. By reducing the rigidity of the compliant blades still 
further, and increasing their unstrained, or "as fabricated”, 
camber, the leading edge loading differences noted in the figure 
could be magnified. However, this optimization process has not been 
carried forward at the time of completing this Final Report. Only 
the potential of the method has been demonstrated. 
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APPENDIX A 



CASCADE AERODYNAMICS 

Initially we wish to derive the steady aerodynamic load on a cascade 
of camberlines in incompressible flow. The usual thin airfoil assump- 
tion is made that the vorticities (representing the jump in tangential 
velocity across the airfoils) are distributed along the chordlines and 
the boundary condition are satisfied along these same straight segments 




Fig. 1 Cascade Geometry 



On the x-axis the velocities induced by the individual airfoils are 
given by* 



U - c 




a 



I 



-lat;- 



1 

2iri 



A-tC 




cL'C. 



The Biot-Savart Law expresses the conjugate velocity u,-iv at a general 
field point. Subsequent applications of the Plemelj formulae supply 
the limiting values of u-iv as the camberline is approached from above 
(+) or below (-) . 
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where ^ -V- A 

pitch, and so on for each airfoil. 

Summing the contribution of all the 
suits in 



and A = S 6 is ths complex 



(infinite number of) airfoils re- 
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Now a well-known result is 



4 CO 



fUs-CD S 



+ LAVlrJ ~ da-’flu ^ 
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so that 






/n=- CO 






where 




I 



Dropping the prime notation this may be put in the standard form 






F(jf'X;T) -t-iCCt-K-,'?) •= y>cAL [p(?-X)3 - (2I>) 

o.r>cL i~ C^) ^ O ^ (y Co) 5= O 
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Since g is purely real, separation of Eq. (2) into real and imaginary 
parts yields 

^ (S 

a ^ ^ X V F ( V><) ^<5 ( ^ 

The boundary condition expresses the fact that the flow direction must 
be tangential to the blade surface (camberline in this case) : 

^ ^ or 4- \/ o(- 

U -F V C<rt3 ^ 

where V and oi. refer to the mean flow and are defined in the figure^^ 
Introduction of Equations (3) and (4) in Equation (5) results in the 
following integral equation 

Solution of Integral Equation 

If the LHS of this equation is presumed defined, then the inversion is 
well known, and is given by 

which may be written 

In the application of the B.C. the term + k g(x) is omitted from W.(x) 
because of a fundamental approximation in thin airfoil theory*. 



Despite any local nonzero slope of the camberline the component «, alone 
changes from one side of the boundary to the other. However, this jump 
in U. is the local magnitude of the vorticity and physically the vortici 
represents the jump in that component of the velocity which is parallel 
to the local boundary. Hence omitting the local vorticity (or strength 
of the vortex sheet) at the point where the B.C. are being applied, the 
remaining components of velocity must sum to a vector parallel to the 
local boundary. Another way of stating this condition 'is to say that 
if all components of the velocity, excepting the local vorticity, sum t< 
a resultant which is parallel to the boundary, the subsequent addition 
of the vorticity component will not affect that parallelism (only the 
magnitude) . 
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where is the so-called single airfoil contribution 

%• |'/5 ri 



(f) 



and the kernel K is given by 






A1 

t->< 



(V 



If a modified vorticity function, f, is defined by 









(1^) 



then the integral equation (7) takes the alternative form 

t 



In order to obtain numerical solutions to this equation the follow- 
ing coordinate transformation, followed by discretization of the 
variables, is introduced. 





X ^ 


5 


A_ 'fjP' 


fi^O 


/f ^ 


n - 


CAd 4>) 






/ 1 \ , 






-f - 




1 ' 


1 - 




i\Xc) 



where ^ ^ and p are positive integers less than or equal to 
or else zero. 



Under this parameterization Eqs . (8) and (9) reduce to 
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4(e)" 4 Q f«'^ 

Ccrj'\^~ C^]9 I J 

and the integral appearing in Eg. (11) reduces to 



If the numerators of the integrals in Eqs. (8' ) and (9* ) are assumed 
to be expanded in cosine series, it may be shown that these Glauert- 
type integrals may be approximated by sums 



/w 




where 






A. 






x±_to,4_ ... 

o , 



'Yl 



V-A^ 






r/- 



f.ir 



/>b 



and S is the Kronecker delta, 
rs 
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In Eq. (13), on the other hand, a direct application of the trape- 
zoidal rule results in 



I nv 



/' 



When these individual results are combined in integral equation 
(11) there results a set of n/1 simultaneous equations, each 
equation corresponding to a particular value of the index ^ . 
These may be expressed in matrix notations as 



^ H 









wherej^Jjis a n/1 square matrix in which ^ identifies the row and 






the column, 



The 



term is 







Both andil/jare column matrices with n/1 rows in which the 
rows are given, respectively, by 





(>7) 



/K. 









(\S) 
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where 





Al 



£(B1 0 



_J 






i(T )('+“- 




(^i) 

^7) 



Knowing the circulation allows the mean flow velocity V and direc 
tiono(to be related to the upstream (or downstream) resultant 
flow. Arranging the cascade along a vertical axis, the various 
velocity vectors take on the aspect shown in the sketch. 




From the geometry of this 
sketch we deduce 



-tan o( = tano^. 
X I 



2 tan ) 






whereas the dimensionless 
circulation for one pitch 
distance along the 
cascade is given by 



sV 






Replacing 
Eg. (25) 



-"6UrO^ 
to Eg. 



■a. in 
(29) 



Eg. (29) 
results 



by Eg. (28) and then eguating 
in 
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Partitioning u into two parts 






“>< — c-ro o' (A 



where 
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The solution of Eq. (15) may be expressed 






d /r 

ji 



i') 



/c ^ 

J- 






where 









At 
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A' 
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Cascade Properties 
The circulation round one blade is 

r = c jVxVx - oVj'^ 



.*7) 

P r t ^ ^ © 



')t 



il~ i£-tix 



4:' 






-5-^ (Q_/i^o<f— 



(oZs') 



Presuming that c><:f may be determined by this means, it is then possible 
to express V in terms of the physically determinate 



This results in the interesting conclusion that a given cascade may 
be solved up to the determination of R and Q without regard to the up- 
stream flow conditions. Subsequently the velocity triangles may be 
deduced for any upstream flow condition by finding from (30) 

and (31) and applying (28). In fact these manipulations may be used 
to derive Merchant's formula 




( 3 / ) 




showing the linear dependence of on o/ • 
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APPENDIX B 



PLATE MECHANICS 



The moment acting on the spring supported plate, see Fig. 2, is 
expressed by first determining the aerodynamic moment about the 
leading edge , 

(m) 



MOM, - c - -k 



This allows the computation of the rearward 
immediately from equilibrium considerations 



spring force 



= MOM. 



(B3) 



Then computation of the a^odynamic lift 

lifi - foU - J (/-*■ 

allows computation of the leading edge spring force 

4c = ~ tC 

Knowing the spring constants k€ and kr the leading edge and 
rearward support deflections may be computed in turn 



y(',n 

^-frAr 



(PA 



where Y is the colum matrix (vector) of elastic displacement 
dimensionless w.r.t. the chord and the indices are given in terms 
of 1 < indix < m s n+1. 

The bending moment distribution may now be computed rearward 
from the ^trailing edge according to . 

{^G) 



Al^fte- X) yji ^ -x) 



—C.<9S 






where U is the unit step function and f , T are dummy variables of 
integration analogous to x, 6. 

In these integrals (Bl, B6 and others) the integrals are 
determined at the discrete mesh points; hence the numerical 
integrations may be evaluated using either a trapezoidal rule, a 
Simpson's rule or other appropriate numerical technique. 

For the approximation to the curvature d^ydx^=cd^Y/dX^ several 
trial models were evaluated. Eventually, in order to perform the 
integration with uniform steps in the angular variable (0 or \|r) the 
following approximation was developed. 
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with V=y/c and Xsx/c=^ (1-cos 0) it is found that 




and expressing the derivatives in terms of control differences, 






where u=jr/m and a common factor h is retained to reduce the size of 
the elements in the D matrix. 

In the application of this result it is found that the first 
and last rows in the (n+l)x(n+l) matrix are not required since the 
bending moment is zero at leading and trailing edges. Recognizing 
that two values of Y(g,l) are known from previous equilibrium 
considerations, two columns may be removed and the useful D matrix 
may be reduced to (n-l)x(n-l) dimensions. In this reduced form the 
solution for the remaining Y terms may be obtained by inverting the 
reduced D matrix and solving for the remaining unknown values of Y 



where the subscript r refers to the reduced versions of these 
matrices. 

Upon reassembling the full Y vector it is then possible to 
estimate the elastic slopes of the camberline using central 
differences* to correct the stagger angle and the camberline slopes 
for the next iteration 



Adding the elastic slopes to the unstrained camberline slopes the 
total slopes are then made available as inputs to the aerodyamics 
subroutine for the next iteration. 



*For leading edge and trailing edge resp. forward and backward 
differences are used. 





^new=^oid+arctan(Y(m, 1) -Y (1, 1) ) 



Y' s z 



Z,^„=tan(arctan + arctan (Y(m, 1) -Y(l, 1) ) 
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APPENDIX B 



PLATE MECHANICS 





Fig. 2 Structural Model with Aerodynamic Load 
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APPENDIX C 



Computer Programs 
test 00. m 
solve. m 

in MatLab language 
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Fig. 3 Computation Flow Chart for solve. m Program (Matlab) 
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Leading Edge Loading 



Rigid and Compliant Blade Comparison 




Fig. 4 Leading-Edge Loading Parameter for 
Rigid and Compliant-Bladed Rotors 
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